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Abstract 

A formulation to include the effects of wall proximity in a second-moment 
closure model that utilizes a tensor representation for the redistribution terms 
in the Reynolds stress equations is presented. The wall-proximity effects are 
modeled through an elliptic relaxation process of the tensor expansion coef- 
ficients that properly accounts for both correlation length and time scales as 
the wall is approached. Direct numerical simulation data and Reynolds stress 
solutions using a full differential approach are compared to the tensor repre- 
sentation approach for the case of fully developed channel flow. 

1. INTRODUCTION 

The theoretical development of higher order closure models, such as Reynolds 
stress models, have primarily been formulated based on high Reynolds number as- 
sumptions. The influence of solid boundaries on these closure models has usually been 
accounted for through either a wall function approach or a modification to the high 
Reynolds number form of the pressure-related correlations and tensor dissipation rate 
and predicated on the near-wall asymptotic behavior of the various velocity second 
moments (So et al. 1991, Hanjalic 1994). 

A broader based attempt to account for the proximity of a solid boundary is 
the elliptic relaxation approach introduced over a decade ago (Durbin 1991) and 
further developed for second-moment closures (Durbin 1993a, Wizman et al. 1996; 
Manceau and Hanjalic 2000, Manceau, Carlson and Gatski 2001). In its two-equa- 
tion form the v 2 — f model has been applied to a variety of flows (e.g., Durbin 1993b, 
1995; Pettersson Reif et al. 1999). The new approach outlined here introduces a ten- 
sor representation for the combined effects of a near-wall velocity-pressure gradient 
correlation and anisotropic dissipation rate that asymptotes to a high Reynolds num- 
ber form away from solid boundaries through an elliptic equation for the polynomial 
expansion coefficients. The development of a generalized methodology for determin- 
ing the polynomial expansion coefficients of representations for the turbulent stress 
anisotropies by (Gatski and Jongen 2000) is extended to an elliptic relaxation proce- 
dure for these expansion coefficients. 

Although the material presented here introduces tensor representations and a ten- 
sor projection methodology into the elliptic relaxation formulation, this work can also 
be viewed as an intermediate step between a fully explicit elliptic relaxation algebraic 
Reynolds stress formulation and the full differential elliptic relaxation Reynolds stress 
formulation. 
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The predictive capabilities of the new model are assessed through comparisons 
with direct numerical simulation channel flow data (Moser et al. 1999). These com- 
parisons include both mean and turbulent flow quantities. 


2. Theoretical Background and Development 

In this section, a mathematical framework is developed for the Reynolds stress 
transport equations and the corresponding elliptic relaxation equation when a tensor 
representation of the redistribution terms is used in the formulation. The method- 
ology introduces a set of elliptic relaxation equations for the polynomial expansion 
coefficients of the chosen representation. The r — f model uses the redistributive 
terms in the elliptic equations, while the r — f3 n model uses the expansion coefficients 
in the elliptic equations. Both models use the Reynolds stress transport equations. 


2.1 Transport Equations 

The transport of the Reynolds stresses r,j (= — uiuj ) is governed by the equation 


Dry, 
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dUi 

= ~Tik^ 

OXk 
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( 1 ) 


where U t is the mean velocity, is the pressure redistribution term, £,j is the tensor 
dissipation rate, and and D'f i are the turbulent transport and viscous diffusion, 
respectively. In the development outlined here, it is best to have (f>ij given by 


(f>ij — m 


dp 

dxj 


dp 2 d pu k 
Uj dx { + 3 dx k 


Su , 


( 2 ) 


so that the trace of the pressure redistribution term is zero. In the application of the 
elliptic relaxation method, it is also necessary to account for the effect of the dissi- 
pation rate anisotropy as the wall is approached. This accounting for the dissipation 
rate anisotropy is accomplished (e.g., Manceau 2000) by a relaxation of the dissipa- 
tion rate anisotropy to its wall value, which is assumed to be equal to the Reynolds 
stress anisotropy. This assumption allows the Reynolds stress transport equation in 
(1 ) to be written as 


£t Vj 
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dUj dUi 
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+ 1 A fij — 
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(3) 


where 

sD fij = 4>ij ~ 2e (dij — bij ) , (4) 

with the Reynolds stress anisotropy 6 i; and dissipation rate anisotropy d-,j defined as 


i _ T ij _ ^ij_ 

ij ~ 2K 3 ’ 


_ $ij 

ij ~ 2e 3 ' 


(5) 


The original scaling of the relaxation function fij was solely through the turbulent 
kinetic energy A'; however. Manceau, Carlson and Gatski (2001) have recently shown 


that an added dissipation rate factor, £, to the scaling (eKfij) eliminates an unwanted 
amplification effect inherent in the original scaling. 

Equation (3) is closed when the model for the turbulent transport T ) \ j is used. In 
previous elliptic relaxation studies that used the Reynolds stress transport equations, 
the viscous diffusion and turbulent transport terms were modeled as 
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where C TJ< — 6 determines the switch to the Kolmogorov time scale (t'/e) 1 ^ so that 
the turbulent time scale will not vanish as the solid boundary is approached. Away 
from the boundary, the composite time scale asymptotes to the inertial scale Kje. 

In the two-dimensional flow considered here, solutions were obtained for the T\\ 
and T22 normal Reynolds stresses and the r 12 shear stress. A transport equation for 
the turbulent kinetic energy was obtained from one-half the trace of Eq. (3) and was 
solved for in lieu of the third normal stress 733, 


,JL(r nn 

Dt + dx t V 'V dx k J 


+ ^V 2 /v, 


( 8 ) 


where V = TikdUi/dxk- The modeled transport equation for the turbulent dissipation 
rate £ needed for closure is given by 


D£_ 

Dt 
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where a e = 1.3, C £ 1 = 1.44, C e 2 = 1.83, with 

C* x = C e i ^1 + a iT^) ’ = 9-09. 


( 10 ) 


Note that this form of the dissipation rate equation (Durbin 1991) has introduced the 
composite time scale into both the production and destruction of dissipation terms. 


2.2 Elliptic Relaxation Methodology: r — f Model 

The rescaled elliptic relaxation equation is driven by the high Reynolds number 
form of the pressure-strain rate correlation II and a contribution from the Reynolds 
stress anisotropy 2£ c bij (away from the wall the dissipation rate is assumed to be 
isotropic dq = 0). This combination results in an elliptic relaxation equation for fij 
given by (cf. Manceau and Hanjalic 2000) 

(l - hi = 7^ K + 2 eMi) - ft, (11) 
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where 


( 12 ) 


and the relaxat ion scales are defined as 
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with Cl = 0.16 and Cl k = 80. Previous implementations of the elliptic relaxation 
procedure (Manceau and Hanjalic 2000) using the Speziale, Sarkar and Gatski (SSG) 
pressure strain rate model (Speziale et al. 1991) used the full nonlinear form. The 
linear form of the SSG model implemented here is given by 

n^j = — (Cj*£ c + C\V^ b{j + A C 2 Sij + A C*3 ^ bikSkj + Sikbkj — -~b nm S nm Sij^j 


—KC A (bikWkj — Wikbkj) (14) 


with Ci = 1.8, Cl = 3.4, C 2 = 0.37, C 3 - 1.25, and C 4 = 0.4. Note that since the 
linear form of the pressure-strain rate model is used here, the value for Cl differs 
from that used previously (Cl = 0.2, see Manceau and Hanjalic 2000) for the form 
of the elliptic relaxation equation given in (11). 

Boundary conditions are needed for the / tJ and are determined, in the vicinity 
of the wall, by the balance of the redistributive term by the viscous diffusion of the 
Reynolds stresses resulting in Table 1. Only the 22- and 12-components of f have 
determinate solutions to the near-wall balance of the stress transport equations. For 
the remaining components f n = fa — —fa/'Z are used as boundary conditions to en- 
sure that fij is traceless (Manceau. Carlson and Gatski 2001). Symmetry conditions 
were applied at the centerline. 

In the current work, one of the goals is to develop a methodology for incorporating 
a tensor representation for the relaxed redistribution function f tJ . Once developed and 
validated this same procedure can be used in conjunction with tensor representations 
for the Reynolds stress anisotropies as well. Such a combination would then yield an 
elliptic relaxation explicit algebraic stress model. The details of the representation 
for the Reynolds stress anisotropy will not be addressed in the current work, but 
deserves further work. As will be discussed in Sec. 2.3, such a representation would 
be consistent with a linear pressure-strain rate model. 


2.3 Representations and Elliptic Relaxation: r — f3 n Model 

Although the elliptic relaxation formulation has already been applied within a full 
differential Reynolds stress model, a question arises about what role tensor represen- 
tations can play within the framework of the elliptic relaxation procedure. 
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Table 1. Boundary Conditions for the fij Tensor 


Component 

Wall 

Centerline 

/ n 

f 22,w 
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/■22 
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— 20// 2 ri2 

i 
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u wy i 

Symmetry 

f 33 

Symmetry 

/ 12 
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The differential elliptic relaxation equat ion for f ZJ can be obtained from the inte- 
gral expression (e.g. Manceau and Hanjalic, 2000) 

' F {j fax') 
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Reynolds stress tensor 


(16) 


and G'n(x. x') is approximated by the free-space Greens function G’n(x.x') = ( 47 rr) _1 
with r = 1 1 x' — x|| . The S ZJ contributions to both the dissipation rate and Reynolds 
stress anisotropies cancel so that the only remaining contributions are the tensor 
dissipation rate and Reynolds stress tensor. The tensor function fij and Tij can be 
represented by polynomial expansions of basis tensors just as the associated Reynolds 
stress anisotropy tensor has been. For such a basis given by T--"\x) (m = 1 , . . . , N), 
the following representations are assumed: 
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fij fa) = £A(x)4 ?) (x) 
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F*(x.x') 


(18) 


ftj(x.x') 

e(x)A’(x) 
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= 5Z 7n(x,x')^ n) (x). 


A tensoi’ scalar product (denoted by [ : ]) between each basis tensor 7’;” * ( x ) and 
the representations given in Eqs. (17) and (18) can be formed, and this leads (using 
matrix notation for convenience) to 


X>(x) [T('>(x):T< m >(x)] = f 


F*(x, x^rT^^x)] G 0 (x,x') 


N . 

= Ey <* 3 x' 7 n(x,x') [T‘^(x):T , ^(x)j6 , n(x,x')(19) 

71=1 ^ 

Since the functional dependency of the indicated scalar product depends solely on x. 
Eq. (19) can be rewritten as 

Ai(x) = f f/ 3 x / 7„(x.x')G'n(x.x / ). (20) 

The modeling of the scalar function 7 n (x. x') follows that established previously for 
the elliptic relaxation approach, that is 

7»(x,x') = 7 „(x\x / )exp • (21) 

where, in general, the 7 n coefficients can have an associated length scale uniquely 
defined by the foi’m given in Eq. (13). 

With this model, Eq. (20) can be rewritten as 

/?n(x)= / cPx' 7 n (x', x 7 ) eXP ^ ^ ■ (22) 

Jn 47r r 

This equation leads directly to the differential counterpart 

(l - L 2 n Y 2 ) fl n (x) = -X*7n(x, x) = /# (x), (23) 

where (3„(x) are the expansion coefficients from the tensor representation of a quasi- 
homogeneous form of f. Since the dissipation rate is assumed to be isotropic, f is 
composed of the quasi-homogeneous form of the pressure-strain rate correlation and 
a contribution due to the Reynolds stress anisotropy. The resultant expression for 
/?n(x) is given by 
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£>£(x) [t< b >(x):T<"*>( X ) 
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e(x)A'(x) 


[(n' l (x) + 2 £c b(x)) :T< m >(x)] 
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(24) 
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where the quasi-homogeneous form of the pressure-strain rate model IT is given by 
n* = £ c (c? - 2 4- C{ ^) b + KC 2 S + KC 3 (bS + Sb - |[b:S]l) 

— KC 4 (bW — Wb) . (25) 

Note that a comparison of Eqs. (14) and (25) shows that the return-to-isotropy term 
proportional to b has been modified. The factor s c now influences the enth'e term and 
the contribution from the Reynolds stress anisotropy 2e c b to the relaxation function 
f is now included in this (slow) term contribution to Il£. 

One of the improvements in the current elliptic relaxation formulation is that 
the scaled relaxation function f defined in Eq. (4) is (9(1) in the log-layer region. 
This scaling negates the adverse influence of the elliptic operator in the log-layer that 
occurred in the original (Durbin 1993a) formulation. In order to retain this benign 
effect in the tensor representation formulation used here, it is necessary to ensure that 
the expansion coefficients 3 n also have this neutral effect. 

Previous representations for the Reynolds stress anisotropy tensor have used basis 
tensors of the form S, SW — WS. and S 2 — [S:S]I/3. In the log-layer, where the 
velocity gradient has a y~ l behavior, this choice of basis tensors would require that 
the corresponding expansion coefficients /?i, /3 2 , and j3 3 have a y, y 2 , and y 2 behavior, 
respectively, in that region to ensure that f behaves as (9(1). Unfortunately, given 
that behavior of the (3 n , the amplification effect would now effect the j3 n and the 
sought-after (9(1) behavior for the f is lost. For the fully developed channel flows of 
interest, this problem can be easily circumvented by using a normalized basis set of 
the form 

T (1) = S‘, T (2) = S*W* - W*S*, T (3) = S* 2 - i, (26) 

where S* = S/jS 2 } 1 / 2 and W* = W/jS 2 } 1 / 2 . This normalization now makes the 
behavior of both the expansion coefficients and basis tensors (9(1) in the log-layer, 
which then precludes any adverse effect of the elliptic operator in the relaxation 
equation (23). 

Boundary conditions for the (3 n expansion coefficients are required. Consistent 
with the boundary conditions for the tensor function fij , the corresponding /3 n bound- 
ary conditions are listed in Table 2 as functions of (see Appendix A for details). 

The equivalence of the elliptic relaxation of the expansion coefficients j3 n given 
by Eq. (23) with the elliptic relaxation of the function fij given by Eq. (11) can be 
readily shown with the current normalized basis. The solution to Eq. (24) is easily 
obtained as 

(si;-. ft, ft) = 777 0 n ‘ :T “'l’ i rT ‘ :T ' 2 'l- [n*:T 121 ]) 

= jr (x/2ft* 12 , - /7‘ ,), 3 (n' c ‘ u + o) . (27) 
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Table 2. Boundary Conditions for j3 n 
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If the tensor representation Eq. (17) is applied to f-j, then the (3 n solution from Eq. 
(27) would yield for the components of /,* 



A comparison of the right-hand side of Eq. (28) with the right-hand side of Eq. (11) 
shows that the two are equivalent. (The reader should recall from the discussion 
following the definition of Il£ in Eq. (25) that the form of the slow term was slightly 
modified from the definition given in Eq. (14). With this change taken into account, 
the exact equivalence Eqs. (11) and (28) holds.) 

3. Results and Discussion 

All flow calculations were carried out on fully developed turbulent channel flows. 
The equations that were solved were scaled in wall units with friction Reynolds 
number Re T based on channel half-height and friction velocity at the wall. A one- 
dimensional finite-difference algorithm described in Appendix B was used for all com- 
putations. 

As shown in Sec. 2.3, the representation methodology that has been developed 
yields an elliptic relaxation formulation equivalent to the elliptic relaxation of the ten- 
sor function While such tensor projection methods have been used in conjunction 
with nonlinear algebraic equations, the application here also validates its use with 
differential operators. 

Figures 1-3 show the predictive accuracy and equivalence of both the — 
fij and Tij - (3 n approaches. The flow field is the fully developed channel flow at 
Re T = 590 (Moser et al. 1999). The figures include both a linear and log scale in 
the wall normal direction. As can be seen from Fig. 1 for the mean velocity, both 
the distribution across the channel and the near-wall asymptotic behavior agree with 
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the direct numerical simulation (DNS) data. Excellent agreement with the DNS data 




Figure 1. Mean velocity distribution across channel at Re T = 590: (a) log-linear scale; 
(b) log-log scale. 

across the channel is also shown for the shear stress profile (Fig. 2); however, the 
asymptotic approach to the wall is greater than the theoretical estimate of 0(y 3 ). 
The discrepancy becomes apparent for values of y < 1. This result is in contrast 
to the predictions for the turbulent kinetic energy shown in Fig. 3. In this case 




Figure 2. Turbulent shear stress distribution across channel at Re T = 590: (a) linear 
scale; (b) log-log scale. 

the near-wall asymptotic behavior is consistent with the DNS results but the overall 
values are slightly lower across the channel than the DNS data. Overall, the predictive 
results for the mean velocity, Reynolds shear stress and turbulent kinetic energy are 
quite exceptional and show that the method can be calibrated to provide excellent 
predictions of this flow field. In actuality, since the models are formally equivalent, 
no changes are required in any of the calibration constants. 

Since a full differential Reynolds stress model is used for the turbulent velocity 
field, it is possible as well as insightful to further examine the component stress 
predictions. Figures 4 and 5 show the r u and t 22 component stresses. Since the 
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Figure 3. Turbulent kinetic energy distribution across channel at Re T = 590: (a) 
linear scale; (b) log-log scale. 

near-wall asymptotic behavior 0(y 2 ) is dominated by the r u (and T33) components, 
it is not surprising to see from Fig. 4b that the near-wall asymptotics closely match 
the DNS results. The 0(y 4 ) behavior that characterizes the DNS results for the T 22 
component (see Fig. 5b) are very closely replicated by the predictions. Fig. 4a shows 
that across the channel predicted results were lower than the DNS results for the T\\ 
component. For the r 2 2 component, however, the predicted peak value was higher 
than the DNS results, but the predicted values were lower over the remainder of the 
channel, as seen in Fig. 5a. 




Figure 4. Reynolds normal stress component ru distribution across channel at Rt T 
= 590: (a) linear scale; (b) log-log scale. 

An interesting assessment of how well the elliptic relaxation formulation models 
the redistribution terms across the channel can be obtained from Eq. (4). The 
quantity eh fij obtained from the explicit representation given in Eq. (18) and the 
elliptically relaxed (3 n from Eq. (23) are plotted in Fig. 6 along with the quantity 
(f>ij— 2 e(dij—bij) obtained from the DNS data. As Fig. 6a shows, the eKf 12 component 
produces the corresponding DNS results very well in the near-wall region and in the 
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Figure 5. Reynolds normal stress component r 2 2 distribution across channel at Re r 
= 590: (a) linear scale; (b) log-log scale. 

outer layer region toward the centerline. Between these two regions, the peak value 
of the computations greatly exceeds that of the DNS. The normal sKfn and e/\/ 2 2 
components show an even poorer prediction of the DNS results. In these cases, only 
the outer layer region is correctly predicted; whereas, over the rest of the channel the 
qualitative and quantitative predictions are generally poor. While the results of this 
a priori validation of the elliptically relaxed function el\ fij are disappointing, it is 
clear that the actual predictions of the fully modeled set of equations are generally 
very good. Thus other modeled terms in the formulation are able to account for any 
discrepancies in the prediction of the redistribution term. 

As Fig. 6 shows, all components of the elliptically relaxed redistribution term 
correctly reproduce the DNS data in the outer layer of the channel flow but differ 
extensively from the DNS data when reproducing the inner layer. Since the elliptic 
operator term (— L 2 V 2 ) is responsible for the deviation of the j3 n from their quasi- 
homogeneous f3„ forms, it is worthwhile to quantify the size of the region across the 
channel that is affected by this term. Figure 7 shows the distribution of —L 2 \ 72 /3 n 
across the channel for the three expansion coefficients ( n = 1,2,3) at three different 
values of Re r . In the inner layer, the wall unit scaling basically collapses the results 
for all values of Re r , with the exception of the component where the results in the 
near- wall region show some dependence on Re T ; this sensitivity to Re T is not found in 
the other components as Figs. 7b and 7c show. The effect of the elliptic operator falls 
to zero at y (wall unit) values around 10 2 . The overshoot in the outer layer shown in 
all the figures is attributed to the asymptotic behavior of the energy dissipation rate 
£ in this region. Both V 2 f3 n and e decrease (L increases); however, the dissipation 
rate e decreases faster ( L increases faster) than the corresponding decrease in V 2 /3 n . 
The variation with Re r in this region is not surprising since the wall unit scaling is 
not the proper scaling for this region. 
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Figure 6. Comparison of predicted redistribution term components with DNS data 
at Re T = 590: (a) 12-component, (b) 11-component, (c) 22-component. For all 
components r 4 j — d n results are eKfij, and DNS results are 
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Figure 7. Effect of elliptic operator across channel at different Re r : (a) /^-coefficient, 
(b) /^-coefficient, (c) /^-coefficient. 
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4. Summary 

A methodology has been developed that introduces a polynomial representation 
for the tensor redistribution function fa. An elliptic relaxation equation, analogous to 
the fa relaxation equation is formulated for the polynomial expansion coefficients /3 n . 
The new prediction method is demonstrated on a fully developed channel flow problem 
and gives similar results to the previous elliptic relaxation method for fa. A formal 
equi valence is established between the elliptic relaxation of the tensor function fa and 
its tensor representation. Although the predictions of the mean velocity and turbulent 
stresses are generally accurate over the channel, an a priori assessment shows that 
the current formulation does not model the redistribution well. Such results are 
enlightening but are not uncommon; the results reflect the fact that in modeled 
closure schemes, a combination of modeled terms combine to yield predictions of 
quantities such as the mean velocity and Reynolds stresses. 

While the theoretical approach developed here does not result in a reduction 
in computational cost, it does introduce a new methodology that is requisite for 
developing elliptic relaxation explicit algebraic stress models. The next step in the 
development of such models will be to introduce representations for the Reynolds 
stress anisotropies and analyze the effects of modeling the turbulent transport and 
viscous diffusion terms consistent with the approximations made in the formulation 
of algebraic stress models. 
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Appendix A fa Boundary Conditions 

The expressions for the fa boundary conditions are derived from the basis tensors 
T/”* used in the representation of fij 


k = E AT 
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(Al) 


Table 1 gives the corresponding boundary conditions for these fij components. The 
boundary condition for (3\ is directly proportional to the /12 boundary condition and 
is given by 


= ^ry = v/2/i^ = 
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The coefficient fa appears in all three expansions of the diagonal terms of fij. If 
fij is traceless, a unique expression for fa at the wall will be obtained. From the 
representation for [ 33 , fa <w can be immediately written as 
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— 30l/ 2 T 2 2 
£ lV(i) 


(A3) 


The representations for f n and can be used to obtain an equivalent expressions 
for the /3s boundary condition 


i3z,u> 


f rp(2) _ r 

J 22jW-*-H 1 w J 11 , tt/ 


7~.(2) 

1 22, w 


7~>(3) 71(2) 7~>(3) 7-1(2) 

^ 22, vu^- 1 l,w * 22,w 


if 22 ,w + fx 


7-1(3) 

^ 33 , w 


— 3/33, ty 


(A4) 


With fa' W known, the representation for either fn or / 22 can be used to obtain fa. w . 
From the / 22 representation, the wall boundary condition on fa is given by 


fa,w 




22, w 


T, 


(3) 


22.W 


fa,w) 


7 - 1 ( 2 ) 

1 22, w 


{j22,w + 3/33 jtt ;T22^) 


j 1 ( 2 ) 

1 22, w 


/ 22 ,U' + 


fz. 


3 ,w 


— 15 u 2 t 2 % 

e wV(l) 


(A5) 
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Appendix B Numerical Solution Methodology 

A one-dimensional finite-difference code was used for all computations. All equa- 
tions were normalized by the bulk viscosity and friction velocity (i.e., wall units). 
The differencing template was node-centered with clustering close to the wall using 
an exponential stretching function. In terms of the scaling used for the channel flow 
calculations, the 500 node grid had the first point at a height of 0.1 wall units. The 
channel Reynolds number Re T determined the channel grid height. The K and e 
equations were implicitly coupled as were the Tjj — fij equations, and for the second 
modeL the — /3 n equations. The variables U, A , £ and r t j were solved in a time- 
dependent mode, while the fij or /?/ equations were not (i.e., At = 0.) All variables 
were updated at each time step. 

In this appendix, the terms with the superscript (r? + 1) denote variables that 
were implicitly solved for and the terms with the superscript (n) were variables used 
explicitly at each iteration. The Reynolds stress equations coupled implicitly with 
either the fij or f3 n equations were solved first, with the momentum, turbulent kinetic 
energy and dissipation rate equations solved second. An updated Tu was used in the 
momentum equation, but the eddy viscosity in the turbulent transport terms of all 
the equations was not updated until after the completion of each time step. Typically 
solutions were re-started from previous turbulent flow calculations. 

The symbol y\ denotes the height of the first node from the wall and e w denotes 
the boundary condition value for £. The discretized form of the governing equations 
are as follows. For the r — f model; 


(n+l) (»). Af 

' xj 1 ij ' 


pW + £ (n) K (n) f {n+l) 


r (n+1) d 

Jl c(») + _ 

A'( n > dy 


Vi 


(n) ' 


1+ T r " 

a K J dy 


(n+l) 


An+l) r(n) 2 ^_r(n+l) _ \ / TT h , o c h 

Aj L dy iJtj ~ £ (n) K (n • 

The boundary conditions were implicitly written for the fij as 
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For the r — /?/ model; 


(n+l) __ (n) 


IJ 


= r^ + At 


P^ + e^KW (jr^ n+1) fy 
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Similarly, the boundary conditions were implicitly written for the /?/ as 
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For U, A and e; 
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The boundary conditions were implicitly written for £ as 
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